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Introduction  and  Problem  Statement 

Micro  air  vehicles  (MAVs)  are  flight  vehicles  with  geometric  dimensions  less 
that  15  cm  and  weight  less  than  100  grams  [1],  The  need  for  such  vehicles  is 
motivated  primarily  by  unmanned  military  and  civilian  missions  that  involve 
closed  spaces  (such  as  buildings)  or  short  distances.  Existing  MAV  designs  may 
be  classified  into  three  categories:  fixed,  rotary,  or  flapping  wing.  While  fixed  and 
rotary  wing  designs  benefit  from  the  existing  technologies,  flapping  wing  designs 
are  fundamentally  different  and  arc  inspired  from  biological  flyers  such  as  insects, 
and  small  birds  and  bats.  In  particular,  flapping  wing  designs  that  operate  at  low 
Reynolds  number  (10‘  <Re<105)  and  low  forward  flight  speed  (7-17  m/s)  have 
received  considerable  attention  due  to  the  exceptional  flight  capabilities  observed 
in  insects  [1-3]. 

A  significant  portion  of  research  on  flapping  wing  MAVs  has  focused  on  the 
mechanisms  that  generate  the  unsteady  aerodynamic  forces.  This  research  [1-7] 
has  identified  the  predominant  force  generating  mechanisms  as  leading  edge 
vortices  (LEVs),  combined  with  wing  rotation,  and  wake  capture  as  a  result  of 
wing-wake  interaction.  Attempts  to  model  the  aerodynamic  loads  have  employed 
two  approaches:  (1)  computational  fluid  dynamics  (CFD)  based  on  solutions  of 
the  Navier-Stokes  (NS)  equations,  and  (2)  approximate  aerodynamic  models 
based  on  potential  flow  solutions.  Descriptions  based  on  CFD  [7-10]  offer  the 
best  resolution  of  the  unsteady  flow  field;  however  these  approaches  arc 
expensive  when  parametric  studies  are  of  interest.  Approximate  models,  on  the 
other  hand,  offer  a  compromise  between  accuracy  and  computational  efficiency; 
therefore  such  approaches  arc  best  suited  for  trend  and  design  type  studies. 
Historically,  simple  steady  and  quasi-steady  models  were  used;  however,  these 
models  consistently  underestimated  forces  generated  by  flapping  wings  [1], 
Unsteady  aerodynamic  formulations  used  for  flapping  wing  problems  can  be 
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classified  as  assumed  (or  prescribed)  wake  and  free  wake  models.  Assumed  wake 
models,  developed  originally  for  aeroelastic  studies  of  fixed  wing  vehicles  [11], 
are  two-dimensional  linear  thin  airfoil  theories  that  are  valid  for  small  plunge 
(bending)  and  pitch  (torsional)  displacements.  Recent  studies  [12,13]  have 
attempted  to  incorporate  the  effect  of  the  LEV  in  assumed  wake  formulations. 
These  approaches  involve  modifying  Theodorsen’s  theory  by  using  a  leading  edge 
suction  analogy  [14]  that  was  originally  derived  for  steady  separated  flow  on  delta 
wings  for  moderate  angles  of  attack  (typically  less  that  40  degrees).  Free  wake 
models  include  unsteady  panel  methods  [10]  and  discrete  vortex  methods  [15- 
17].  These  approaches,  which  account  for  evolution  of  the  wake,  provide  a 
reasonable  approximation  to  the  development  of  the  unsteady  wake  during  a 
flapping  cycle.  A  two-dimensional  panel  method  was  used  in  Ref.  [10],  and 
reasonably  good  agreement  was  found  with  CFD  simulations  for  the  range  of 
parameters  considered.  However,  LEVs  were  not  accounted  for  in  this  study.  Two 
dimensional  discrete  vortex  formulations  that  account  for  flow  separation  have 
been  considered  in  Refs.  [15-17].  The  model  developed  in  Ref.  [15],  which 
accounts  for  separation  close  to  the  leading  edge,  compared  well  with 
experimental  data  for  airfoils  in  steady  flow.  In  this  approach  the  chordwise 
location  of  the  separation  point,  which  may  be  obtained  using  independent 
computations  or  experiments,  is  explicitly  incorporated  into  the  formulation. 
Comprehensive  treatment  of  an  unsteady  aerodynamic  model  based  on  a  discrete 
vortex  method  that  is  applicable  to  insect-like  flapping  wings  in  hover  is 
presented  in  Refs.  [16,17],  The  model  was  used  to  simulate  rigid  wings,  and  for 
the  cases  considered,  compared  well  with  experimental  data  on  flapping  wings. 

The  importance  of  wing  flexibility  in  enhancing  the  performance  of  flapping 
wings  has  been  mentioned  in  a  number  of  studies;  however,  only  a  limited 
number  of  studies  [13,18,19]  have  attempted  to  address  this  issue  in  a  systematic 
manner.  In  the  earliest  study,  where  a  linear  finite  element  model  of  the  wing  was 
coupled  with  an  unsteady  panel  method  [18],  the  wing  model  was  based  on 
experimentally  obtained  geometric  and  inertial  data  on  Moth  wings.  More  recent 
computational  studies  have  considered  wing  models  based  on  Eulcr-Bemoulli 
beam  theory  coupled  with  CFD  [20],  or  membranes  reinforced  by  metal  or 
composite  frames  [13]  coupled  with  Theodorsen’s  theory  modified  using 
Polhamus  leading  edge  suction  analogy.  These  studies  [13,18,19],  which 
considered  wing  flexibility  in  a  linear  manner,  concluded  that  wing  flexibility  had 
a  favorable  effect  on  lift,  and  that  the  effect  of  flexibility  cannot  be  ignored  when 
computing  aerodynamic  loads. 

Thus,  the  overall  objective  of  the  research  carried  out  in  the  framework  of 
the  current  grant  was  to  determine  the  effect  of  flexibility  on  the  performance  of 
anisotropic  flapping  wings  in  hover  and  forward  flight.  The  specific  objectives 
were  as  follows: 

(1)  Develop  a  nonlinear  acroclastie  model  to  study  flexible  anisotropic 
flapping  wings. 

(2)  Conduct  systematic  validation  studies  to  develop  confidence  in  the  model. 

(3)  Explore  the  effect  of  wing  flexibility  on  aeroelastic  response  and  unsteady 
loads  generated  by  flapping  wings. 
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Nonlinear  Aeroelastic  Model 

The  aeroelastic  model  is  obtained  by  coupling  MARC,  a  commercially  available 
nonlinear  finite  element  code,  with  the  approximate  aerodynamic  model  that  was 
developed  in  Refs.  [16,17];  a  description  of  this  model  is  presented  in  Ref.  [20]. 
Structural  Dynamic  Model 

The  structural  dynamic  model  of  a  bio-inspired  wing,  depicted  in  Fig.  1, 
{note:  all  figures  are  provided  at  the  end  of  the  report )  consists  of  beam  and 
shell  elements  capable  of  undergoing  large  rigid  body  rotations  as  well  as 
moderate-to-large  flexible  deformation.  The  elements  have  compatible  degrees  of 
freedom  so  that  beam  stiffened  shell  structures  may  be  constructed  [21], 
Furthermore,  the  elements  support  a  variety  of  constitutive  laws  so  that  isotropic 
as  well  as  anisotropic  wings  can  be  modeled.  Wing  kinematics  consists  of  large 
amplitude,  prescribed,  time  dependent  rigid  body  rotations  imposed  at  the  base 
(root)  of  the  wing.  In  MARC,  rotations  may  be  prescribed  either  using  a  time 
dependent  rotation  vector  [21]  or  as  displacements  at  two  or  more  nodes. 
Aerodynamic  Model 

The  approximate  aerodynamic  model  is  based  on  potential  flow  and  uses  a 
circulation/vorticity  approach  to  compute  the  aerodynamic  loads.  The  formulation 
accounts  for  leading  edge  separation  and  subsequent  vortex  formation, 
incorporates  the  effect  of  wing  thickness  and  camber,  and  includes  a  free  wake 
model.  The  model  is  two-dimensional  and  is  applied  in  a  strip  theory  manner. 
Furthermore,  each  airfoil  cross-section  is  assumed  to  interact  only  with  its  own 
shed  wake.  At  each  time  step,  vorticity  shed  into  the  wake  is  computed  by 
enforcing  a  Kutta  condition  at  the  trailing  edge  and  a  stagnation  condition  at  the 
leading  edge.  The  evolution  of  the  wake  is  governed  by  the  Rott-Birkhoff 
equation,  which  is  derived  from  the  Biot-Savart  law  for  two-dimensional  flow. 
The  unsteady  aerodynamic  loads  are  computed  using  the  vortex  impulse  method 
and  the  unsteady  Bernoulli  equation.  It  is  important  to  note  that  the  formulation 
[16,17]  was  originally  developed  for  rigid  wings  for  the  case  of  hover;  its 
implementation  for  wings  that  have  both  span-wise  and  chord-wise  flexibility  is 
described  in  Ref.  [20], 

Coupled^  Fluid-Structure  Model 

The  aeroelastic  response  is  obtained  using  an  updated  Lagrangian  (UL) 
method  [22]  wherein  the  equations  of  motion  at  each  time  step  are  formulated  by 
using  the  equilibrium  configuration  computed  at  the  previous  time  step  as  the 
reference.  A  schematic  description  of  the  approach  is  shown  in  Fig.  3.  The 
coupled  fluid-structure  problem  is  modeled  in  MARC  via  user-defined  load 
subroutines  [21]  whereby  time  dependent  pressure  loads  that  are  computed  from 
the  wing  motion  arc  imposed  on  the  structure.  The  user  subroutine  is  called  from 
the  main  program  for  each  Newton-Raphson  iteration  of  the  UL  method  thereby 
ensuring  convergence  of  the  structural  displacements  and  aerodynamic  loads 
within  each  time  step.  Schematic  of  the  aeroelastic  solution  in  MARC  is  depicted 
in  Fig.  4.  Finally,  the  resulting  equations  of  motion  are  integrated  using  a  suitable 
numerical  scheme. 
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Results  and  Discussion 


Validation  o£ the  structural  dynamic  model 

Important  considerations  in  the  structural  modeling  of  MAV  wings  are  (1) 
implementation  of  prescribed,  large  amplitude,  rigid  body  rotations  that  are 
representative  of  wing  kinematics,  and  (2)  accurate  capture  of  the  dynamie  or 
centrifugal  stiffening  effect  in  flexible  structures  undergoing  prescribed  rotary  or 
flapping  motion.  Note  that  all  the  results  that  follow  were  obtained  by  integrating 
the  equations  of  motion  using  a  single  step  Houbolt  scheme  [23], 

Implementation  of  large  amplitude  rigid  body  rotations  in  MARC  was 
examined  by  imposing  kinematics  of  a  dragonfly  wing  on  a  rigid  rectangular  plate 
shown  in  Fig.  2.  The  results  [20],  shown  in  Fig.  5,  indicated  that  wing  kinematies 
can  be  accurately  implemented  as  displacements,  but  produee  significant  error 
when  implemented  as  a  time  dependent  rotation  veetor. 

The  effect  of  centrifugal  stiffening  was  examined  by  considering  the  spin-up 
motion  of  a  flexible  plate;  details  of  the  test  ease  are  given  in  Ref  [24],  The 
results,  shown  in  Fig.  6,  indicate  that  this  effect  is  accurately  modeled  in  MARC. 
Validation  of_ the  aerodynamic  model 

The  aerodynamic  model  was  validated  for  (1)  eases  of  attached  flow  over  the 
airfoil,  and  (2)  eases  where  flow  separation  from  the  leading  edge  was  observed. 

The  force  coefficients  for  a  NACA0012  airfoil  undergoing  prescribed  plunge 
motion  are  shown  in  Fig.  7.  Figure  9  shows  a  comparison  of  lift  coefficients 
obtained  for  a  rigid  flat  plate  airfoil  undergoing  prescribed  plunge  motion  in  near 
hover  conditions.  Parameters:  Chord  =  1.0  m,  plunge  frequency  =  0.064  Hz, 
Reynolds  number  =  1000,  plunge  amplitude  =  0.5m.  The  CFD  results  were 
obtained  by  implementing  laminar  NS  equations  in  CFD++,  a  commercially 
available  CFD  eode,  using  a  grid  that  had  approximately  1 50,000  cells  with  240 
points  on  the  airfoil.  Results  using  the  approximate  model  were  obtained  by 
assuming  flow  separation  at  the  leading  edge.  This  result  indicates  that  the 
approximate  model  shows  reasonable  agreement  with  CFD  for  the  ease 
considered. 

Aeroelastic  response  results 

Preliminary  aeroelastic  results  that  examined  the  effeet  of  flexibility  on  force 
generation  by  wings  undergoing  prescribed  motion  have  been  presented  in  Ref. 
[20].  The  calculations  were  performed  by  assuming  leading  edge  separation  for  a 
zero  free  stream  velocity,  thereby  simulating  conditions  of  hover.  The  results, 
which  were  obtained  for  wings  that  had  different  spanwise  stiffness  [20],  were 
consistent  with  the  finding  of  previous  studies  and  indicated  that  flexibility  had  a 
comparatively  small  but  favorable  impact  on  force  generation  and  that  the  effeet 
of  aerodynamic  loads  on  wing  deformation  was  small  compared  to  the  effeet  of 
inertia  loads.  Sample  results  of  lift  and  chord-normalized  tip  displacements  are 
shown  in  Fig.  10. 

Concluding  Remarks  and  Accomplishments 

The  research  performed  during  the  period  of  this  grant  resulted  in  the 
development  and  comprehensive  testing  of  a  nonlinear  aeroelastic  model  that  is 
suitable  for  the  analysis  of  flexible  anisotropic  MAV  wings  for  the  ease  of  hover; 
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this  work  [20]  was  presented  at  the  49th  AIAA/ASME/ASCE/AHS/ASC 
Structures,  Structural  Dynamics  and  Materials  Conference  held  in  April  2008  in 
Schaumburg,  1L. 

Currently,  the  model  is  being  extended  in  the  frame  work  of  the  multi¬ 
university  research  initiative  (MURI)  to  include  effects  of  forward  flight,  finite 
span  and  tip  vortices.  Furthermore,  a  separation  and  re-attachment  criterion  will 
also  be  included  so  that  the  aeroelastic  model  can  be  used  for  both  hover  and 
forward  flight  for  a  range  of  angles  of  attack. 
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Figure  1:  Structural  dynamic  model  of  a  bio-  Figure  2:  Rectangular  plate  used  for 

inspired  MAV  wing  preliminary  studies 


Equations  of  motion  (EOM)  in 
integral  form  (Principle  of  Virtual  Work) 
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Figure  3:  Schematic  of  the  updated  Figure  4:  Schematic  of  the  aeroelastic  solution 

Lagrangian  method  in  MARC 
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Figure  5:  Implementation  of  wing  kinematics:  Figure  6:  Effect  of  centrifugal  stiffening.  Results 

solid  line  -  MATLAB  computation,  dashed  line  -  from  Jinyang  and  Jiazhen  presented  in  Ref.  [24]. 
rotation  vector  in  MARC,  circles  -displacements 
in  MARC 
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Figure  7:  Lift  and  moment  coefficients  for  a  NACA0012  airfoil  undergoing  prescribed  plunge 

motion.  UPM  results  were  obtained  from  Ref  [10].  k  and  h  denote  reduced  frequency  and  chord- 
normalized  plunge  amplitude  respectively. 
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Figure  9:  Lift  coefficients  for  a  rigid  flat  plate  airfoil  undergoing  prescribed  plunge  motion  in  near 
hover  conditions. 


Figure  10:  Lift  and  chord-normalized  tip  displacement  of  wings  undergoing  prescribed  plunge 
motion  for  the  case  of  hover  [20]  Black  -  rigid  wing,  Red  and  Blue  -  flexible  w  ings;  Blue  line 
corresponds  to  the  most  flexible  configuration  tested. 
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